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A  DOMAIN  DECOMPOSITION  ALGORITHM  BASED  ON 
A  CHANGE  TO  A  HIERARCHICAL  BASIS 

BARRY  F.  SMITH  ♦  AND  OLOF  B.  WIDLUNDf 

Abstract.  Over  the  last  five  years,  several  new  fast  algorithms  have  been  developed  for  the 
solution  of  elliptic  finite  element  problems  with  many  degrees  of  freedom.  Among  them  are  Yserentant's 
hierarchical  basis  multigrid  method  and  a  number  of  domain  decomposition  algorithms  for  problems 
divided  into  many  subproblems.  The  condition  number  of  the  relevant  operators,  for  the  best  of  these 
preconditioned  conjugate  gradient  methods,  grows  only  slowly  with  the  size  of  the  problems;  the  growth 
typically  is  quadratic  in  the  logarithm  of  the  number  of  degrees  of  freedom  associated  with  a  subregion 
or  an  element  of  the  coarsest  mesh. 

Important  components  of  many  domain  decomposition  preconditioners  are  subproblems  related  to 
the  sets  of  variables  on  the  interfaces  between  neighboring  subregions.  In  this  paper,  we  demonstrate 
that,  for  problems  in  two  dimensions,  a  very  simple  change  of  basis  for  these  one-dimensional  problems 
leads  to  a  preconditioner  which  is  as  effective  as  those  previously  considered.  In  our  proof,  we  use  only 
tools  of  linear  algebra  and  a  result  of  Yserentant's.  Our  numerical  experiments  confirm  that  the  new 
algorithm  is  better  conditioned  than  the  original  hierarchical  basis  multigrid  method. 

Key  Words,  domain  decomposition,  elliptic  equations,  finite  elements,  iterative  substructuring, 
hierarchical  bases,  preconditioned  conjugate  gradients,  Schur  complement 
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1.  Introduction.  We  consider  second  order,  self  adjoint,  uniformly  elliptic  dif- 
ferential equations  on  a  two  dimensional  polygonal  domciin  fi.  The  problems  are  solved 
numericcdly  by  using  continuous,  piecewise  linear  finite  elements.  The  domain  is  first 
subdivided  into  nonoverlapping,  triangular  subregions  fi,,  also  called  substructures, 
and  these  are  further  triangulated  into  elements.  H  denotes  the  diameter  of  a  typical 
substructure  and  h  the  diameter  of  its  elements. 

We  develop  a  domain  decomposition  algorithm  similar  to  those  considered  by 
Bj0rstad  and  W^idlund  [3], [4],  Bramble,  Pasciaic,  and  Schatz  [5], [6],  Dryja  and  Widlund 
[13]  and  Widlund  [19].  When  using  these  methods  the  variables  interior  to  individual 
substructures  are  first  eliminated.  The  resulting  reduced  system,  the  Schur  comple- 
ment, therefore  only  involves  the  variables  associated  with  F,  the  set  of  edges  and 
vertices  of  the  substructures.  This  system  is  then  solved  by  a  preconditioned  conju- 
gate gradient  method,  where  the  preconditioner  is  constructed  from  certain  problems 
associated  with  the  interfaces  r,j  =  dQ,iC\dClj  between  the  substructures,  and  a  global 
coarse  problem  associated  with  the  vertices  of  the  substructures.  We  note  that  it  is 
shown  in  Dryja  and  Widlund  [13],  that  such  a  preconditioner  naturally  can  be  viewed 
in  terms  of  a  splitting;  cf.  Varga  [18].  In  the  splitting,  the  couplings  between  the 
groups  of  variables,  associated  with  individual  edges  of  the  substructures,  are  elimi- 
nated. In  [13],  we  also  explain  how  results  and  algorithms  for  the  two  substructure 
case  can  be  used  to  construct  and  analyze  problems  on  many  substructures. 

Various  preconditioners  have  been  proposed  for  the  subproblems  associated  with 
the  edges  r,_,.  For  each  F,j,  this  is  essentially  a  two  subregion  problem,  and  we  can 
therefore  take  advantage  of  a  number  of  results  obtained  in  early  work  on  domain 
decomposition  algorithms.  Alr;eady  in  1980,  Dryja  [11],  see  also  [12],  introduced  an 
effective  preconditioner  J,  which  is  the  square  root  of  a  discrete,  one  dimensional 
Laplacian  on  F,j.  The  same  preconditioner  was  also  discussed  in  Bj0rstad  and  Wid- 
lund [4]  and  Bramble,  Pasciak  and  Schatz  [5].  Other  preconditioners  for  these  two 
subregion  subproblems,  such  as  the  Neumann-Dirchlet  algorithm,  were  considered  by 
Bj0rstad  and  Widlund  [4],  Bramble,  Pasciak  and  Schatz  [6],  Chan  and  Resasco  [8], 
Chan  and  Keyes  [7],  Dihn,  Glowinski,  Periaux  [10]  and  Golub  and  Mayers  [14].  A 
number  of  the  resulting  algorithms  for  the  many  substructure  case  are  known  to  be  al- 
most optimal  in  the  sense  that  the  condition  number  is  bounded  by  C(l  +  \og(H/h)y. 
For  a  more  complete  discussion,  see  Bj0rstad  and  Widlund  [4]  and  Widlund  [19]  for 
the  two  subregion  and  many  subregion  cases,  respectively. 

An  alternative  almost  optimal  algorithm,  which  uses  a  hierarchical  basis,  has  been 
introduced  by  Yserentant  [21].  His  bound  on  the  condition  number  is  of  the  same  form. 
When  the  standard  finite  element  nodal  basis  is  replaced  by  a  hierarchical  basis,  the 
transformed  matrix  becomes  much  better  conditioned.  The  preconditioner  is  then 
obtained  by  discarding  the  off  diagonal  blocks  and  by  replacing  all  but  one  of  the 
diagonal  blocks  by  diagonal  matrices. 

In  this  paper,  we  consider  a  hybrid  method  demonstrating  that  a  successful  and 
simple  preconditioner  can  be  obtained  by  changing  the  bases  of  the  spaces  associated 
with  individual  edges  r,j.  Our  proof  uses  only  tools  of  linear  algebra  and  Yserentant 's 
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main  result.  We  show  that  the  new  method  has  a  smaller  condition  number  than 
Yserentant's  original  method;  thus  it  grows  no  faster  than  C(l  +\og(H/h)Y ,  a  result 
confirmed  in  our  numerical  experiments. 

Our  work  has  been  inspired  by  recent  work  of  Babuska,  Craig,  Mandel  and 
Pitkaranta  [1],  Babuska,  Griebel,  and  Pitkaranta  [2]  and  Mandel  [17],  [15],  [16]  where 
efficient  preconditioners  for  the  p-version  of  the  finite  element  method,  are  devel- 
oped by  using  hieraxchiceil  basis  functions  and  partial  orthogonalization  of  the  basis 
functions.  Similarly,  our  algorithm  for  the  h-version  involves  a  change  of  basis.  It  is 
extremely  easy  to  carry  out,  and  it  results  in  a  much  better  conditioned  linear  system. 

2.   Yserentant's  Algorithm  and  Iterative  Substructuring  Methods.  In 

this  section  we  provide  the  necessary  background  to  make  it  possible  to  define  the 
new  algorithm  and  to  put  it  in  some  perspective. 

We  consider  a  second  order,  self  adjoint,  uniformly  elliptic,  bilinear  form  an(u,  v) 
on  Q  and,  for  simplicity,  impose  a  homogeneous  Dirichlet  condition  on  dO.; 

a{u,v)^{f,v),      \/veH'^{n),       ueHli^). 

For  the  two  levels  of  triangulations  into  substructures  f2,  and  elements  introduced 
earlier,  we  assume  shape  regularity  and  that  they  satisfy  the  usual  rules  of  finite 
element  triangulations;  see  e.g.  Ciaxlet  [9].  V^{^)  and  V^{Q)  are  the  spaces  of 
continuous,  piecewise  linear  functions,  on  the  two  triangulations,  which  vanish  on  the 
boundary  dQ. 

The  discrete  problem  is  then  of  the  form 

(1)  a{uH,  V,)  =  (/,  vh),    Vu,  e  v\n),     UH  e  v\n). 

2.1.  The  Hierarchical  Basis  Method.  The  hierarchical  basis  method  pro- 
vides a  general  purpose  preconditioner  for  second  order  elliptic  problems  in  the  plane; 
see  Yserentant  [21], [20].  The  algorithm  is  given  in  terms  of  a  set  of  spaces  V^' ,i  = 
0,  •  •  •  ,7,  which  are  successive  refinements  by  a  factor  of  two  of  V^°  =  V^ .  V^'  is  the 
set  of  piecewise  linear  finite  element  functions  after  i  levels  of  refinement  from  the 
original  coarse  triangulation  with  V^'  =  V^.  V^  is  a  direct  sum  of  subspaces 

where  V^  =  V^'  \  V'^'-^ .  In  other  words,  V^''  is  the  set  of  piecewise  linear  functions  in 
V''* ,  which  are  zero  at  the  nodal  points  of  the  triangles  of  all  coarser  triangulations,  see 
Figure  1  for  a  one-dimensional  case.  For  the  spaces  V^'',  we  choose  a  basis  of  standard 
nodal  functions  of  V^'  associated  with  the  new  nodes.  The  resulting  basis  for  the 
entire  space  V^  is  much  closer  to  being  orthogonal  in  the  H^  sense,  than  the  standard 
nodal  basis  functions,  and  the  stiffness  matrix  is  therefore  much  better  conditioned. 

Yserentant's  preconditioner  is  a  block  diagonal  matrix  in  the  new  basis.  The  first 
block  is  defined  by  the  finite  element  model  for  the  subspace  V^  and  the  others  are 
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Fig.  1.  Hierarchical  Basis  Functions  jn  One  Dimension 

diagonal.  In  matrix  notation  the  resulting  system,  which  is  solved  by  a  conjugate 
gradient  method,  is  of  the  form 

or  equivalently, 

D-'H'^KHx  =g. 

Here  K  is  the  original  stiffness  matrix,  H  represents  the  transformation  from  the 
hierarchical  to  the  nodal  basis,  and  D  is  block  diagonal  and  obtained  from  H^KH, 
as  described  above.  In  an  implementation,  there  is  no  need  to  represent  the  stiffness 
matrix  explicitly  in  the  new  basis.  Instead,  we  perform  the  basis  change  on  the  vectors 
as  needed. 

In  [21],  Yserentant  develops  the  algorithms  needed  for  performing  the  basis  change 
from  hierarchical  basis  to  nodal  basis  and  back  and  demonstrates  that  each  requires 
fewer  than  2n  additions  and  n  divisions  by  2,  where  n  is  the  dimension  of  the  finite 
element  space.  The  following  algorithm  is  valid  for  both  one  and  two  dimensions. 

Algorithm  to  form  x  <—  Hx 

for  fc  =  1  to  number  of  levels 
for  i  on  level  k 

X,  =  X,  +  (x/i,  +X/2,)/2 

next  i 
next  k 

The  integer  arrays  II  and  72  contain  pointers  to  the  two  neighbors  of  i  which  are  on 
the  next  coarser  level.  We  can  regard  the  algorithm  as  defining  a  factored  form  of 
the  matrix  H .  The  nodal  to  hierarchical  transformation  is  similar.  We  note  that  the 
coarse  mesh  needs  not  be  uniform;  see  Figure  1.  If  the  refinements  are  not  uniform, 
then  the  weights  in  the  algorithms  have  to  be  adjusted.  It  is  also  possible  to  continue 
the  refinement  only  in  selected  subregions. 

2.2.  Iterative  Substructuring  Methods.  Iterative  substructuring  algorithms 
use  a  different  splitting  of  the  space  V'^  into  A'^  +  1  subspaces; 


V' 


K 


haT 
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e  Vo\nN) 


For  each  substructure  fi,-,  we  thus  have  a  subspace  ^(^(Q,)  =  V'TlH^iQ,).  The 
elements  of  V/fa^m  ^^e  piecewise,  discrete  harmonic  functions,  i.e.  they  are  orthogonal, 
in  the  sense  of  the  bilinear  form  a(-,  •),  to  all  the  other  subspaces.  It  is  easy  to  show 
that  an  element  of  V^^^^  is  uniquely  determined  by  its  values  on  T  =  \J^Q.^. 

In  a  first  step  of  many  substructuring  algorithms,  the  variables  interior  to  the 
Cli  are  eliminated.  We  partition  the  vector  x  =  {xi,xb)  and  the  stiffness  matrix 
K  accordingly.  The  system  that  remains  to  be  solved  is,  after  a  block  Gaussian 
elimination  step, 

(2)  Schur{K)xB  —  9- 

Here  Schur[K)  is  a  Schur  complement  defined  by 

Schur{K)  =  Kbb  -  KJqKJiKib. 

A  particular  iterative  substructuring  method  is  defined  by  the  choice  of  a  pre- 
conditioner  for  equation  (2).  Finally,  when  accurate  enough  values  on  F  have  been 
computed,  the  values  elsewhere  are  determined  by  solving  A'^  separate  Dirichlet  prob- 
lems on  the  individual  substructures.  We  note  that  it  is  not  necessary  to  compute  the 
elements  of  Schnr(K)  since,  in  the  conjugate  gradient  iteration,  this  matrix  is  needed 
only  in  terms  of  matrix- vector  products.  Such  a  product  can  be  found  at  the  expense 
of  solving  one  problem  on  each  of  the  substructures. 

A  number  of  preconditioners  can  be  described  as  follows:  We  first  carry  out 
a  partial  change  of  basis,  associating  the  standard  basis  functions  of  V^  with  the 
vertices  of  the  substructures.  In  the  new  basis,  we  represent  the  Schur  complement 
as 


EE     Sev 
EV    ^vv 


Here  Svv  denotes  the  part  of  the  Schur  complement  associated  with  the  vertices  of 
the  substructures  and  See  the  part  associated  with  the  edges  between  substructures, 
etc. 

The  preconditioner  for  this  system  is  given  by 

See      0 
0       Svv 

where  Svv  is  the  coarse  mesh  finite  element  problem  and  See  is  a  block  diagonal 
matrix.  Each  of  its  blocks  is  associated  with  the  variables  of  a  single  edge  F.j.  The 
operator  J,  mentioned  before,  can  be  used  for  this  purpose;  for  other  examples  of 
such  algorithms,  see  the  references  given  in  the  third  paragraph  of  Section  1. 

3.  The  Hybrid  Algorithm.  We  could  combine  the  two  main  ideas  of  Section 
2  as  follows:  We  first  represent  the  stiffness  matrix  in  the  hierarchical  basis  and  then 
eliminate  the  interior  variables  of  all  the  substructures.   We  proceed  by  solving  the 
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remaining  Schur  complement  system  approximately  without  further  preconditioning. 
Finally  we  use  the  resulting  values  as  boundary  data  for  the  local  problems  on  the 
individual  substructures. 

In  the  new  algorithm,  we  proceed  differently,  but  as  we  will  see,  we  will  obtain  the 
same  approximate  solution  without  the  considerable  expense  of  converting  the  stiffness 
matrix  into  the  hierarchical  basis.  In  our  algorithm  we  work  with  the  standard  nodcil 
basis  while  eliminating  the  interior  variables,  only  changing  to  the  hierarchical  basis 
on  r,  the  set  of  interfaces  and  vertices.  The  resulting  linear  system  is  similar  to  that 
of  Section  2.2 

Dg^gHlsSchur{K)HBBiB  =  9 

It  is  important  to  note  that  we  do  not  use  any  further  preconditioning  of  the  variables 
associated  with  the  edges  r,j. 

This  method  offers  several  possible  advantages  over  standard  hierarchical  basis. 
The  conjugate  gradient  iteration  is  carried  out  over  a  much  smaller  set  of  unknowns 
and  we  will  show  that  the  condition  number  is  smaller.  The  solution  of  the  subprob- 
lems  is  easily  parallelizable  since  they  are  independent.  The  hierarchical  bases  method 
in  its  original  form  appears  to  offer  less  opportunity  for  this  trivial  type  of  paralleliza- 
tion.  The  change  of  basis  required  in  each  iteration  step  now  consists  of  completely 
independent  one  dimensional  problems  instead  of  a  two  dimensional  problem.  The 
basic  observation  is  that  the  values  at  a  node  on  r,_,  can  be  computed  using  only  the 
coefficients  for  the  hierarchical  basis  functions  related  to  that  edge. 

We  now  prove  the  almost  optimality  of  our  algorithm  using  two  simple  lemmas 
and  Yserentant's  result. 

Lemma  1.  Let  G  represent  a  change  of  basis  which  leaves  the  space  of  variables 
on  r  invariant.  Then  the  Schur  complement  associated  with  this  set  of  unknowns  is 
independent  of  the  choice  of  bases  for  VQ(Q.k)- 

Proof  Let  xj  be  the  vector  of  unknowns  associated  with  VQ{Qk),  Vt,  and  xb  be 
those  associated  with  T.  The  most  general  basis  transformation  considered  here  is  of 
the  form 

Gn     Gib 
0      Gbb 

In  the  new  basis,  the  stiffness  matrix  is 

K=l    ^n       ^     \  (  ^^''f     ^''B  \  (  Gu     Gib 


Gib    Gbb  J  \  -^^/b    ^bb  J  \    ^      Gbb 
A  straightforward  calculation  shows  that  its  Schur  complement  satisfies 


Schur(K)  =  GggSchur(K)G 


BE 


The  following  result  follows  easily  by  a  Rayleigh  quotient  argument. 
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Refinement  levels 12  3  4  5 6 7_ 

#  of  unknowns  on  T        3  7  15  31  63  127  255 


New  method  1.68     2.66      3.82       5.18       6.75        8.52        10.50 


No  preconditioning       3.05     6.88     14.20     28.63     57.37     114.79     230.49 


Table  1 
Condition  Number  for  Two  Subdomain  Case 


Lemma  2.  Let  K  he  symmetric,  positive  definite.  Then,  the  condition  numbers 
of  K  and  its  Schur  com,plem,ent  satisfy 

K{Schur{K))  <  k(A'). 

Our  main  result  is  given  in 

Theorem  1.  The  condition  number  of  the  hybrid  algorithm,  introduced  in  this 
section,   is  bounded  by  that  of  Yserentant's  method.     Thus,   it  is  bounded  by  C(l  + 

\og{H/h)r. 

Proof  We  use  Lemma  1  twice  and  the  fact  that  D  is  block  diagonal  to  obtain 

Schur{D-''^H'^KHD-^l^)    =    D^'i'' Schur{H^  KH)DbT 

=    Ds'J'HlBSchur{K)HBBDsT  ■ 

By  using  Lemma  2,  we  obtain 

K{D:BTHlBSchur{K)HBBDsT)    =    K{Schur{D-'/^H'^KHD-'^')) 

<  k{d-^i^h'^khd-'I% 

which  is  bounded  by  C(l  +  log(i///i))^  see  Yserentant  (Thm.  4.1)  [21]. 

I 

4.  Numerical  Experiments.  In  a  first  set  of  experiments,  we  consider  the  do- 
main Q  =  Qi  U  f^2  where  Cii  and  f22  are  unit  squares  aligned  along  an  edge  F  = 
fii  n  ^2-  We  use  the  standard  regular  mesh  and  the  usual  five  point  discretization  for 
the  Laplacian.  The  results  are  listed  in  Table  1. 

Remark.  Our  experiments  show  that  the  condition  number  grows  faster  than  (1  + 
log  H/h)  for  the  two  subdomain  case.  We  note  that  for  a  number  of  preconditioners 
the  condition  number  remains  bounded  in  this  case.  This  is  true  for  the  preconditioner 
J  if  we  solve  a  Dirichlet  problem,  cf.  [4],  but  not  for  a  Neumann  problem.  Yet,  our 
method  and  that  based  on  the  J  operator  both  have  condition  numbers  which  grow 
like  (1  +  log  if /h)^  in  the  many  subdomain  case. 

In  a  second  set  of  experiments,  we  consider  the  case  of  many  substructures.  The 
unit  square  Q  is  subdivided  uniformly  into  4,  16,  64  or  256  square  subdomains  and 
the  same  model  problem  is  solved  using  uniform  meshes.  We  compare  our  results 
with  a  set  of  experiments  reported  in  Yserentant  [21].  His  coarse  problem  has  only  a 
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Levels 

3 

4 

5 

6 

7 

8 

9 

jj^  of  unknowns  in  Q. 

72 

152 

31^ 

63^ 

1272 

2552 

511^ 

Yserentant's  method 

New  method 

No  preconditioning 


4  Subdomains 
10.59     19.53     31.85     47.14 
3.66       6.13       9.37      13.41 
9.77      21.50     44.97     91.98 


65.38     86.51     110.49 


#  of  unknowns  on  F 

13          29          61 

125 

16  Subdomains 

New  method 

4.85       7.94 

11.80 

16.44 

No  preconditioning 

35          75 

155 

316 

#  of  unknowns  on  T 

81         177 

369 

753 

64  Subdomains 

New  method 

5.29 

8.52 

12.54 

17.32 

No  preconditioning 

137 

290 

599 

1217 

#  of  unknowns  on  T 

385 

833 

1729 

3521 

256  Subdomains 

New  method 

5.39 

8.70 

12.78 

17.60 

No  preconditioning 

546 

1152 

2372 

4766 

#  of  unknowns  on  F 


1665      3585      7425      15105 


Table  2 
Condition  Numbers  for  the  Many  Subdomam  Case 


single  node  and  it  is  therefore  directly  comparable  only  with  our  experiments  on  four 
substructures.  The  results  are  given  in  Table  2. 

The  coarse  problem  and  the  problems  associated  with  the  edges,  which  together 
make  up  the  preconditioner,  are  independent.  We  can  therefore  scale  the  contribution 
of  the  coarse  model  by  a  scalar  factor  a  selecting  the  value  of  the  parameter  for  which 
the  convergence  is  fastest.  In  our  numerical  experiments,  we  have  found  that  for  our 
model  problem  a  ss  3.6  is  the  best  for  a  wide  range  of  refinements.  In  the  case  of  4 
subdomains  with  a  coarse  problem  consisting  of  a  single  point,  we  have  found  that  the 
best  Q  varied  considerably  with  the  mesh,  ranging  from  5.4  to  7.89.  We  note  that  the 
condition  number  grows  quadratically  in  the  logarithmic  factor  for  any  fixed  a  >  0. 
Our  numerical  results  are  reported  for  optimal  choices  of  <y. 

In  other  experiments,  which  will  be  reported  elsewhere,  we  have  found  that  there 
is  very  little  difference  in  the  performance  of  the  new  and  several  previously  known 
domain  decomposition  methods.  We  also  note  that  a  version  of  our  algorithm  has 
been  tested  successfully  for  a  membrane  model  of  plane  linear  elasticity  by  Anders 
Hvidsten  and  Petter  Bj0rstad  of  the  University  of  Bergen,  Norway. 
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